Evaluating land consolidation decision for smallholder farmers in Kenya
Hannah Kamau, Cory Whitney
1 Introduction
This document is a supplementary resource and contains the description of the land consolidation model used to simulate outcomes of the decision to consolidate land in Kenya by smallholder farmers. Data used in this model, was collected in March 2024 after co-creating a conceptual model with the stakeholders.
2 Install the necessary packages
Must have R packages are decisionSupport(Luedeling et al. 2024),
tidyverse(Wickham 2023), and
ggplot2 (Wickham 2016).
#Packages citation tied to a .bib file
knitr::write_bib(c(.packages(),
'decisionSupport',
'ggplot2',
'tidyverse',
'magrittr',
'kableExtra',
'DT',
'grid',
'scales'), 'bib/packages.bib')3 Transcribing conceptual model to a land consolidation function
Conceptual model co-created with stakeholders reflected the costs, benefits, and risks that could arise when shifting from the small and highly diversified farms to a model where farmers consolidate their land under the management of someone else. Varied benefits and costs can be envisioned from the current system of small farms and from a model where farmers consolidate land. For example,
Costs that a farmer could incur in a land consolidated option
- Planning costs (legal fee, knowledge acquisition)
- Loss of natural assets e.g., trees
- Food costs not saved
- Household upkeep cost prior to the first payment
Benefits that a farmer could get in a land consolidated option
- Lease income
- Saved production costs
- Alternative income from on or off farm activities
- Better childhood
- Social cohesion
- Saved medical bills related to land use accidents/diseases
Costs incurred by a farmer who remains in the current system of small farms, henceforth referred to as “do nothing” option
- Planning costs (Whether the farmer decides to change or not, (s)he will have incurred the cost of planning)
- Medical bills related to land use accidents/diseases
- Production costs
Benefits that accrue in the do nothing option
- Land security (calculated by the value of trees)
- Annual crop yield
- Saved food costs
The above options are also faced by risks. Risks facing farmers who shift include social risks as a result of freed time while risks associated with crops will face farmers in the do nothing option. NB: Shifting to a land consolidation model, is for a period of 25 years. During this time frame, the farmer can longer farm on the piece of land consolidated, but can receive lease income commensurate to the size of the land.
Below is a series of equations that define the above costs, benefits,
and risks for both options in the function farmer. The
calculations help us find out what land consolidation truly renders by
calculating the difference between the options. This difference is the
DECISION that can inform on viability of land consolidation for the
farmers. To calculate the decision, the below equation takes a wide
range of inputs, which means, uncertainty is unavoidable.
set.seed(254)
farmer <- function(x, varnames)
{
# A) do nothing benefits ####
# We use maize as the reference crop
# ex-ante risks to farmers crops
natural_hazard
pest_disease_risk
farmer_inadequate_funds_risk
late_planting_risk
health_risk
prop_maize_yield_lost_hazard <- chance_event(natural_hazard,
yield_loss_drought/100, 0, gen_CV, n=n_years,
one_draw = FALSE)
prop_maize_yield_lost_disease <- chance_event(pest_disease_risk,
yield_loss_disease/100, 0, gen_CV,n=n_years,
one_draw = FALSE)
prop_maize_yield_lost_input_constraint <-
chance_event(farmer_inadequate_funds_risk, yield_loss_due_to_input_constraint,
0, gen_CV, n=n_years, one_draw = FALSE)
prop_maize_yield_lost_management <-
chance_event(min(late_planting_risk, health_risk),
yield_loss_due_to_management, 0, gen_CV, n=n_years,
one_draw = FALSE)
#effects of all the above risks
farm_risks <- sapply(c(prop_maize_yield_lost_hazard +
prop_maize_yield_lost_disease + prop_maize_yield_lost_input_constraint
+ prop_maize_yield_lost_management),function(x) min(x,1))
farmer_yield_t_ha <- vv(maize_yield_t_ha, gen_CV, n_years) * no_of_seasons
adjusted_farmer_yield_t_ha <- farmer_yield_t_ha * (1 - farm_risks)
crop_benefit <-
adjusted_farmer_yield_t_ha *
ha_per_hh * 1000 * # conversion to kg/ha
vv(maize_price_kes_kg, gen_CV, n_years)
# natural capital
adjusted_value_of_assets <-vv(value_of_farm_assets, gen_CV, n_years,relative_trend = inflation_rate)
# Food cost saved from food accessible from the farm
annual_saved_food_cost <- vv(saved_food_cost_pm, gen_CV, n_years,relative_trend = inflation_rate) * 12
status_quo_benefit <- crop_benefit + adjusted_value_of_assets + annual_saved_food_cost
# B) do nothing costs ####
# medical Bills from farm related stresses
medical_bills <- prop_hhincome_spent_on_hospital/100 * vv(hh_income_pa, gen_CV,n_years)
# production costs for the farm
farming_costs <- vv(production_costs_saved_acre, gen_CV, n_years,relative_trend = inflation_rate)*
ha_per_hh * ha_acre_conversion * no_of_seasons
# planning cost involves legal fee and the cost of knowledge acquisition
farmer_plan_cost <- planning_cost
status_quo_cost <- medical_bills + farming_costs + farmer_plan_cost
# C) Farmer costs with land consolidation ####
farmer_one_time_cost <- cost_of_disruption_kes + # Damages on physical assets
hhupkeep_prior_to_first_payment # HH income needed to sustain the hh prior 1st payment
# planning cost involves legal fee and the cost of knowledge acquisition
farmer_plan_cost <- planning_cost
farmer_recurring_cost <- adjusted_value_of_assets + # loss
annual_saved_food_cost # Food unavailable from farm directly
land_consolidation_cost <- farmer_one_time_cost + farmer_plan_cost + farmer_recurring_cost
# D) Farmer Benefits with land consolidation ####
# farmer's annual compensation
lessor_fee <-vv(compensation_income_pm_acre, gen_CV, n_years)* 2 * ha_per_hh * ha_acre_conversion
#off farm employments and wages
off_farm_income <-vv(off_farm_income_kes_pm, gen_CV, n_years) *12
#income saved: production costs for farming
production_costs_saved <- farming_costs
# hospital bills saved (proxy for long term health) from farm related stresses
medical_bills_saved <- medical_bills
# social cohesion benefit
# social cohesion is expressed as a factor of free time as a result of the the intervention. With the intervention, free time is associated with social cohesion through community participation and involvement. We quantify it by the cost of labor per hour. First, risks that threaten social cohesion are anti-social behaviors e.g crime, drug use, theft and domestic conflict - They take up time that would otherwise be used for community cooperation or activities
social_time_risk <- max(vice_risk,domesticconflict_risk)
adjusted_social_time <- vv(social_time_hr_day, gen_CV, n_years) * (1-social_time_risk)
social_cohesion <- adjusted_social_time * vv(labour_cost_kes_hr, gen_CV, n_years,
relative_trend = inflation_rate)
# Better childhood benefit = long term benefit is education quantified by child's contribution to family farm labor weekly and the value of hired labor
better_childhood <- child_farm_time_hr_week * children_per_hh * n_weeks_child_farm_time *
vv(labour_cost_kes_hr, gen_CV, n_years, relative_trend = inflation_rate )
land_consolidation_benefit <-lessor_fee + off_farm_income + production_costs_saved +
sale_of_hh_items_not_needed + medical_bills_saved +
social_cohesion + better_childhood
# E) Calculate net benefit ####
status_quo_netbenefit <- status_quo_benefit - status_quo_cost
status_quo_netbenefit_usd <- status_quo_netbenefit/currency_change
land_consolidation_netbenefit <- land_consolidation_benefit - land_consolidation_cost
land_consolidation_netbenefit_usd <- land_consolidation_netbenefit/currency_change
# Calculate NPV categorized costs and benefits ####
food_cost_saved_npv <- discount(annual_saved_food_cost, discount_rate,calculate_NPV = TRUE)
natural_assets_npv <- discount(adjusted_value_of_assets, discount_rate,calculate_NPV = T)
lease_npv <- discount(lessor_fee, discount_rate, calculate_NPV = T)
prdn_costs_npv <- discount(farming_costs, discount_rate, calculate_NPV = T)
alt_income_npv <- discount(off_farm_income, discount_rate, calculate_NPV = T)
better_childhood_npv <- discount(better_childhood, discount_rate,calculate_NPV = T)
social_npv <- discount(social_cohesion, discount_rate, calculate_NPV = T)
medical_npv <- discount(medical_bills_saved, discount_rate, calculate_NPV = T)
crop_npv <- discount(crop_benefit, discount_rate, calculate_NPV = T)
# NPV ####
NPV_land_consolidation_netbenefit <- discount(land_consolidation_netbenefit, discount_rate, calculate_NPV = TRUE)
NPV_land_consolidation_netbenefit_usd <- discount(land_consolidation_netbenefit_usd, discount_rate, calculate_NPV = T)
NPV_status_quo_netbenefit <- discount(status_quo_netbenefit, discount_rate, calculate_NPV = TRUE)
NPV_status_quo_netbenefit_usd <- discount(status_quo_netbenefit_usd, discount_rate, calculate_NPV = T)
# Benefit cost ratio ####
# without intervention
npv_status_quo_benefit <- discount(status_quo_benefit, discount_rate, calculate_NPV = T)
npv_status_quo_cost <- discount(status_quo_cost, discount_rate, calculate_NPV = T)
bcr_status_quo <-npv_status_quo_benefit/npv_status_quo_cost
# with intervention
npv_land_consolidation_benefit <- discount(land_consolidation_benefit, discount_rate, calculate_NPV = T)
npv_land_consolidation_cost <- discount(land_consolidation_cost, discount_rate, calculate_NPV = T)
bcr_land_consolidation <- npv_land_consolidation_benefit/npv_land_consolidation_cost
return(list(benefit_lc_usd = NPV_land_consolidation_netbenefit_usd,
benefit_lc_kes = NPV_land_consolidation_netbenefit,
benefit_status_quo_usd = NPV_status_quo_netbenefit_usd,
benefit_status_quo_kes = NPV_status_quo_netbenefit,
net_benefit_lc_kes = NPV_land_consolidation_netbenefit - NPV_status_quo_netbenefit,
net_benefit_lc_usd = NPV_land_consolidation_netbenefit_usd - NPV_status_quo_netbenefit_usd,
BCR_status_quo = bcr_status_quo,
BCR_lc = bcr_land_consolidation,
Food_money_saved = food_cost_saved_npv,
Natural_assets = natural_assets_npv,
Planning = farmer_plan_cost,
Disruption_cost = farmer_one_time_cost,
Lease = lease_npv,
Production_costs = prdn_costs_npv,
Alternative_income = alt_income_npv,
Childhood = better_childhood_npv,
Social = social_npv,
Medical = medical_npv,
Annual_yield = crop_npv,
Cashflow_decision_do = land_consolidation_netbenefit - status_quo_netbenefit))
}4 Monte Carlo simulation
Now we run a Monte Carlo simulation of the model function above using
the inputs form the input table
land_consolidation_input.csv (shows the variables used in
the model function and their input given as ranges - lower and upper
bound, and the type of distribution).
mcSimulation_results <- decisionSupport::mcSimulation(
estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep =";"),
model_function = farmer,
numberOfModelRuns = 1e3,
functionSyntax = "plainNames")5 Results of the simulation
5.1 Distribution of the outcomes
y_output<- as.data.frame(mcSimulation_results$y)
x_output <- as.data.frame(mcSimulation_results$x)
mc_output <-cbind(y_output, x_output)
mc_output <- mc_output %>% mutate(NPV_per_ha = benefit_lc_kes/ha_per_hh,
No_NPV_per_ha = benefit_status_quo_kes/ha_per_hh,
NPV_per_ha_decision = net_benefit_lc_kes/ha_per_hh )
fig1a <- mc_output %>%
select(NPV_per_ha, No_NPV_per_ha) %>%
pivot_longer(., cols = c("NPV_per_ha", "No_NPV_per_ha"),
names_to = "variables", values_to = "project_outcome") %>%
group_by(variables) %>%
mutate(lower_bound = quantile(project_outcome, 0.05, na.rm = T),
upper_bound = quantile(project_outcome, 0.95, na.rm = T)) %>%
filter(project_outcome >= lower_bound & project_outcome <= upper_bound) %>%
ungroup() %>%
mutate(variables = factor(variables,
levels=c('NPV_per_ha', 'No_NPV_per_ha'))) %>%
ggplot(aes(x = project_outcome/1e6, colour = variables, fill = variables))+
geom_density(alpha = 0.3)+
theme_bw()+
xlab("NPV per hectare in million KES") + ylab("Density")+
ggtitle("a) Overall outcomes") +
scale_colour_manual(name = "",
labels = c("Land consolidation","Status quo"),
values = c("#1a80bb", "deeppink"))+
scale_fill_manual(name = "",
labels = c("Land consolidation","Status quo"),
values = c("#1a80bb", "deeppink"))+
theme(legend.position = c(0.98, 0.98),
legend.justification = c("right", "top"),
legend.background = element_rect(fill = "transparent"))
# Decision land consolidate
percentages <- mc_output %>%
select(net_benefit_lc_kes) %>%
summarise(
negative = sum(net_benefit_lc_kes < 0, na.rm = TRUE) / n() * 100,
positive = sum(net_benefit_lc_kes >= 0, na.rm = TRUE) / n() * 100)
negative_percent <- round(percentages$negative, 2)
positive_percent <- round(percentages$positive, 2)
custom_colors <- c("negative" = "#1a80bb", "positive" = "deeppink") #pink
fig1b <- mc_output %>%
mutate(outcome_category = ifelse(net_benefit_lc_kes >= 0, "positive","negative")) %>%
select(net_benefit_lc_kes, outcome_category) %>%
pivot_longer(., cols ="net_benefit_lc_kes",
names_to = "variables", values_to = "project_outcome") %>%
group_by(variables) %>%
mutate(lower_bound = quantile(project_outcome, 0.05, na.rm = T),
upper_bound = quantile(project_outcome, 0.95, na.rm = T)) %>%
filter(project_outcome >= lower_bound & project_outcome <= upper_bound) %>%
ungroup() %>%
ggplot(aes(x=project_outcome/1e6, fill = outcome_category)) +
geom_histogram(breaks=seq(-50,50, by=2)) +
scale_fill_manual(values = custom_colors) +
guides(fill = FALSE)+
labs(x=" NPV in Million KES", y = "Frequency") +
annotate(geom = "text",label = paste0(negative_percent, " %"), x = -34, y = 37,color = "black") +
annotate(geom = "text", label = paste0(positive_percent, " %"), x = 20, y = 37,color = "black") +
theme_bw() + ggtitle("b) Decision outcome for land consolidation") +
theme(plot.margin = margin(0, 0, 90, 0),
axis.title.x = element_text(margin = margin(t = 5))) +
coord_cartesian(clip = "off") + # Allow annotations outside plot area
annotation_custom(grob = linesGrob(arrow = arrow(length = unit(0.25, "cm"), type = "open"), #arrow facing right
gp = gpar(col = "black")), xmin = 5, xmax = 15, ymin = -15, ymax = -15) +
annotation_custom(grob = textGrob("Land consolidation \npreferable",
x = unit(0.8, "npc"), y = unit(-0.1, "npc"),
just = "left", gp = gpar(fontsize = 10)), xmin = 22, xmax = 22, ymin = -15, ymax = -15 ) +
annotation_custom(grob = linesGrob(x = unit(c(1, 0), "npc"), # Reverse direction for arrow
arrow = arrow(length = unit(0.25, "cm"), type = "open"),
gp = gpar(col = "black")), xmin = -5, xmax = -15, ymin = -15, ymax = -15) +
annotation_custom( grob = textGrob("Do nothing \npreferable",
x = unit(0.8, "npc"), y = unit(-0.1, "npc"),
just = "left", gp = gpar(fontsize = 10)), xmin = -50, xmax = -50, ymin = -15, ymax = -15 )
fig1a + fig1b
Figure 1. Distribution of the outcomes given as net present value and expressed per hectare of land. (a) the distribution of land consolidation and do nothing options (b) distribution of the DECISION for land consolidation. The outcomes are based on 1,000 model runs generated through Monte Carlo simulation.
# ggsave("figures/fig_1_MC_outcome_dist.png", width = 9, height = 5.5)5.2 Variable of Importance and Expected Value of Perfect Information
pls_result <- plsr.mcSimulation( object = mcSimulation_results,
resultName = names(mcSimulation_results$y[7]), ncomp = 1)
fig2a <- decisionSupport::plot_pls(pls_result,
input_table = read.csv("Input_tables/land_consolidation_input.csv", sep=";"),
threshold = 1,
pos_color = "#1a80bb", #blue
neg_color = "#ea801c") + #orange
scale_y_discrete(labels = function(y) str_wrap(y, width = 10)) +
ggtitle("a) Variable in Importance") +
theme(axis.text = element_text(size = 9),
axis.title = element_text(size = 12),
legend.text = element_text(size = 10))
mcSimulation_table <- data.frame(mcSimulation_results$x, mcSimulation_results$y[1:6])
evpi <- decisionSupport::multi_EVPI(mc = mcSimulation_table,
first_out_var = "benefit_lc_usd")## [1] "Processing 6 output variables. This can take some time."
## [1] "Output variable 1 (benefit_lc_usd) completed."
## [1] "Output variable 2 (benefit_lc_kes) completed."
## [1] "Output variable 3 (benefit_status_quo_usd) completed."
## [1] "Output variable 4 (benefit_status_quo_kes) completed."
## [1] "Output variable 5 (net_benefit_lc_kes) completed."
## [1] "Output variable 6 (net_benefit_lc_usd) completed."fig2b<- plot_evpi(evpi, decision_vars = "net_benefit_lc_kes",
input_table = read.csv("Input_tables/land_consolidation_input.csv", sep = ";"),
bar_color = "#1a80bb",
x_axis_name = "Expected Value of Perfect Information KES",
new_names = "") +
scale_y_discrete(labels = function(y) str_wrap(y, width = 20)) +
ggtitle("b) Value of Information") +
theme(axis.text = element_text(size = 8),
axis.title = element_text(size = 10))
fig2a + fig2b
Figure 2. Sensitivity and value of information analyses of decision land consolidation. a) Important variables identified by PLS regresson analysis and expressed as VIP scores to indicate model sensitivity. b) Key uncertainties identified by expected value of perfect information metric. The outcomes are based on 1,000 model runs generated through Monte Carlo simulation.
ggsave("figures/fig_2_MC_VIP_EVPI.png", width = 9, height = 5)6 Monte Carlo simulation with scenarios
Based on two most critical variables as shown in Figure 2b, maize and land prices, we created scenarios to report on how the decision changes with varying prices of maize and lease prices. To define the updates of maize and land price variables, we relied on different sources of data to create the bounds of the prices of each. For example, for maize prices, we looked at the daily maize price submissions on the Ministry of Agriculture and Livestock Department on https://amis.co.ke/site/market. Here, we took data published between 2005 to 2024 for Murang’a County.
folder_path <- "Input_tables/maize_prices/"
csv_files <- list.files(path = folder_path, pattern = "*.csv", full.names = TRUE)
list_of_dfs <- lapply(csv_files, function(file) {
read.csv(file, sep = ";") # read the csv files into data frame.
})
# Optionally, add a column with the file name to track the source of each data frame
list_of_dfs <- lapply(seq_along(list_of_dfs), function(i) {
df <- list_of_dfs[[i]]
df <- df %>%
mutate(Supply.Volume = as.numeric(Supply.Volume)) # Convert to numeric (or as.character() if preferred)
df$source_file <- basename(csv_files[i]) # Add file name as a new column to track the source
return(df) })
combined_df <- bind_rows(list_of_dfs)
combined_df <- combined_df %>%
mutate(
Wholesale = gsub("/Kg", "", Wholesale), # Remove "/Kg"
Wholesale = gsub(",", "", Wholesale), # Remove commas (if they exist)
Wholesale = trimws(Wholesale), # Remove any leading/trailing whitespaces
Wholesale = na_if(Wholesale, " -"), # Convert " - " to NA
Wholesale = as.numeric(Wholesale), # Convert to numeric
Retail = gsub("/Kg", "", Retail), # Same process for Retail
Retail = gsub(",", "", Retail),
Retail = trimws(Retail),
Retail = na_if(Retail, " -"),
Retail = as.numeric(Retail),
Date = as.Date(Date)) # Change the date to be date
country <- combined_df %>%
select(Wholesale, Retail) %>%
pivot_longer(., cols = everything(), names_to = "market_type", values_to = "prices") %>%
ggplot(aes(x=factor(market_type), y= prices))+
geom_boxplot(width = 0.1, color = "forestgreen") +
xlab("Market type") + ylab("Price in KES/Kg") +
theme_minimal() +ggtitle("a) Prices of dry maize in the country")
muranga <-combined_df %>% filter(County == "Muranga") %>%
select(Wholesale, Retail) %>%
pivot_longer(., cols = everything(), names_to = "market_type", values_to = "prices") %>%
ggplot(aes(x=factor(market_type), y= prices))+
geom_boxplot(width = 0.1, color = "forestgreen") +
xlab("Market type") + ylab("Price in KES/Kg") +
theme_minimal() + ggtitle("b) Dry maize prices in Murang'a county")
country+muranga
Figure.S1. Distribution of dry white maize prices in Kenya (a) and in Murang’a country (b).
Within the decisionSupport is a function called
scenario_mc, that allows one to create scenarios using
specific variables. one has to create a separate .csv file with specific
variables to be modified. Based on uncertain variables -maize and land
prices- we created market prices scenario file
market_scenarios.csv. The file is as below
Table 1. Specific model input variables that need to be updated to generate the market scenarios. These variables are chosen based on the analysis of the value of information and the results of EVPI are as shown in Figure 2b. These variables are the most critical variables for the decision to consolidate land.
market_scenario_param <- read.csv("Input_tables/market_scenarios.csv", sep = ";")
market_scenario_param %>%
datatable(extensions = "Buttons",
options = list(dom = "Blfrtip",
buttons=c('copy', 'csv', 'excel'),
scrollX = TRUE,
autoWidth = TRUE),
rownames=FALSE)Running the farmer model function with the updated .csv
file which we refer to as Monte Carlo simulation with scenarios
land_consolidation_scenarios <-scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = read.csv("Input_tables/market_scenarios.csv", sep =";"),
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")6.1 Varying the price of maize per kg
NB: When the decision outcomes are positive (i.e., have positive values), land consolidation is preferred and when the outcomes are negative, do nothing option is preferable.
scen_2 <- scenarios[,c(1,2,4)]
scen2_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_2,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
fig3a <-ggplot(scen2_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision (Mio KES)") + theme_bw()+
ggtitle("a) KES 5 - 86.67")
scen_3 <- scenarios[,c(1,2,5)]
scen3_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_3,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
fig3b <-ggplot(scen3_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision (Mio KES)") + theme_bw()+
ggtitle("b) KES 86.68 - 168.33")
scen_4 <- scenarios[,c(1,2,6)]
scen4_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_4,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
fig3c <-ggplot(scen4_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision (Mio KES)") + theme_bw()+ ggtitle("c) KES 168.34 - 250")
fig3a + fig3b + fig3c
Figure S2. The decision outcome of land consolidation when maize price is varied. The decision outcome is based on 1,000 model runs of land consolidation model function. Mio = Million.
6.2 Varying the price of land per acre per season
scen_5 <- scenarios[,c(1,2,7)]
scen5_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_5,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen5 <-ggplot(scen5_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("a) KES 10K -340k")
scen_6 <- scenarios[,c(1,2,8)]
scen6_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_6,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen6 <-ggplot(scen6_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("b) KES 340K -670k")
scen_7 <- scenarios[,c(1,2,9)]
scen7_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_7,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen7 <-ggplot(scen7_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("c) KES 670K - IM")
scen5 + scen6 + scen7
Figure S3. The decision outcome of land consolidation when land price is varied. Outcome based on 1000 model runs of land consolidation model function. Mio = Million.
6.3 Combinations of maize and land prices
scen_8 <- scenarios[,c(1,2,10)]
scen8_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_8,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen8 <-ggplot(scen8_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("a) low MP & LP")
scen_9 <- scenarios[,c(1,2,11)]
scen9_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_9,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen9 <-ggplot(scen9_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("b) low MP & medium LP")
scen_10 <- scenarios[,c(1,2,12)]
scen10_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_10,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen10 <-ggplot(scen10_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("c) low MP & high LP")
scen_11 <- scenarios[,c(1,2,13)]
scen11_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_11,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen11 <-ggplot(scen11_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("d) medium MP & low LP")
scen_12 <- scenarios[,c(1,2,14)]
scen12_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_12,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen12 <-ggplot(scen12_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("e) medium MP & LP")
scen_13 <- scenarios[,c(1,2,15)]
scen13_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_13,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen13 <-ggplot(scen13_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("f) medium MP & high LP")
scen_14 <- scenarios[,c(1,2,16)]
scen14_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_14,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen14 <-ggplot(scen14_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("g) high MP & low LP")
scen_15 <- scenarios[,c(1,2,17)]
scen15_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_15,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen15 <-ggplot(scen15_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("h) high MP & medium LP")
scen_16 <- scenarios[,c(1,2,18)]
scen16_consolidation <- scenario_mc(base_estimate = estimate_read_csv("Input_tables/land_consolidation_input.csv", sep=";"),
scenarios = scen_16,
model_function = farmer,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
scen16 <-ggplot(scen16_consolidation$y, aes(net_benefit_lc_kes/1000000))+
geom_histogram(aes(y=..density..), color="black", fill = "white")+
geom_density(alpha =.3, fill = "#FF8")+
ylab("Relative probability") +
xlab("Decision outcome (Mio KES)") + theme_bw()+
ggtitle("i) high MP & LP")
(scen8 + scen9 + scen10)/
(scen11 + scen12 + scen13)/
(scen14 +scen15 +scen16)
Figure S4. The decision outcome of land consolidation when maize and land prices are varied. These outcomes are based on 1000 model runs of land consolidation model function. Mio = Million.
7 Effect change of the critical variables
Here, we compare the net benefit of land consolidation of different scenarios to the baseline net benefit of land consolidation. Table 2 shows the summary statistics of the net benefit for land consolidation, which can also be referred to as the decision of each scenario.
land_consolidation_scenarios$y[,"net_benefit_lc_per_ha"] <- land_consolidation_scenarios$y$net_benefit_lc_kes/land_consolidation_scenarios$x$ha_per_hh #standardise to per ha
land_consolidation_scenarios$y[,"net_benefit_lc_per_ha_MKes"] <-land_consolidation_scenarios$y$net_benefit_lc_per_ha/1000000 # in millions for simplicity
mcSimulation_results$y[,"net_benefit_lc_per_ha"] <- mcSimulation_results$y$net_benefit_lc_kes/mcSimulation_results$x$ha_per_hh
mcSimulation_results$y[,"net_benefit_lc_per_ha_MKes"] <- mcSimulation_results$y$net_benefit_lc_per_ha/1000000
baseline_lc <- mcSimulation_results$y$net_benefit_lc_per_ha_MKes
low_mp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_2"]
med_mp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_3"]
high_mp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_4"]
low_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_5"]
med_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_6"]
high_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_7"]
low_mp_lp <-land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_8"]
low_mp_med_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_9"]
low_mp_high_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_10"]
med_mp_low_lp <-land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_11"]
med_mp_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_12"]
med_mp_high_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_13"]
high_mp_low_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_14"]
high_mp_med_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_15"]
high_mp_lp <- land_consolidation_scenarios$y$net_benefit_lc_per_ha_MKes[land_consolidation_scenarios$x$Scenario == "Scen_16"]
scenarios_summary <- rbind(data.frame("Baseline" = round(quantile(baseline_lc, c(0,0.05, 0.5, 0.95, 1)),2),
"Low maize prices"=round(quantile(low_mp, c(0,0.05,0.5,0.95,1)),2),
"Medium maize prices"=round(quantile(med_mp,c(0, 0.05,0.5,0.95,1)),2),
"High maize prices"=round(quantile(high_mp, c(0,0.05, 0.5, 0.95,1)),2),
"Low lease prices"=round(quantile(low_lp, c(0,0.05, 0.5,0.95,1)),2),
"Medium lease prices"=round(quantile(med_lp, c(0,0.05,0.5,0.95,1)),2),
"High lease prices"=round(quantile(high_lp, c(0,0.05,0.5,0.95,1)),2),
"Low maize and lease prices"=round(quantile(low_mp_lp, c(0,0.05,0.5,0.95,1)),2),
"Low MP and medium LP"=round(quantile(low_mp_med_lp, c(0,0.05,0.5,0.95,1)),2),
"Low MP and high LP"= round(quantile(low_mp_high_lp, c(0,0.05,0.5,0.95,1)),2),
"Medium MP and low LP"=round(quantile(med_mp_low_lp, c(0,0.05,0.5,0.95,1)),2),
"Medium MP and LP"=round(quantile(med_mp_lp, c(0,0.05,0.5,0.95,1)),2),
"Medium Mp and high LP"=round(quantile(med_mp_high_lp, c(0,0.05,0.5,0.95,1)),2),
"High MP and low LP"=round(quantile(high_mp_low_lp, c(0,0.05,0.5,0.95,1)),2),
"High MP and medium LP"=round(quantile(high_mp_med_lp, c(0,0.05,0.5,0.95,1)),2),
"High MP and LP"=round(quantile(high_mp_lp, c(0,0.05,0.5,0.95,1)),2)))
scenarios_summary <- t(scenarios_summary)
scenarios_summary_df <- as.data.frame(scenarios_summary)
quantile_labels <- c("Minimum", "5% quantile", "Median", "95% quantile", "Maximum")
colnames(scenarios_summary_df) <- quantile_labels
rownames(scenarios_summary_df) <- gsub("\\.", " ", rownames(scenarios_summary_df)) # Replace periods with spaces in variable names
scenarios_summary_df <- cbind("Market price scenarios" = rownames(scenarios_summary_df), scenarios_summary_df, row.names=NULL)
kable(scenarios_summary_df, format = 'html',
caption = "<p style='color:black'><i><b>Table 2.</b> Summary statistics of the net benefit for land consolidation of different market price scenarios. Based on 1000 model runs. Net benefit given in million KES</i></p>",
align = c("l", rep("c", ncol(scenarios_summary_df) - 1))) %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover")) %>%
add_header_above(c("Market Price Scenarios" = 1, "Quantiles" = ncol(scenarios_summary_df) -1), underline = TRUE) %>%
row_spec(nrow(scenarios_summary_df), hline_after = TRUE)| Market price scenarios | Minimum | 5% quantile | Median | 95% quantile | Maximum |
|---|---|---|---|---|---|
| Baseline | -32382.28 | -92.15 | -6.82 | 10.45 | 320.14 |
| Low maize prices | -25901.25 | -93.37 | 0.49 | 14.61 | 7484.82 |
| Medium maize prices | -1398.79 | -81.89 | -6.00 | 9.62 | 327.92 |
| High maize prices | -30201.19 | -103.07 | -12.59 | 2.20 | 227.62 |
| Low lease prices | -2921.71 | -78.35 | -4.70 | 12.72 | 124.47 |
| Medium lease prices | -51621.57 | -47.49 | 15.16 | 36.41 | 70.76 |
| High lease prices | -9188.39 | -36.20 | 32.69 | 54.95 | 236.76 |
| Low maize and lease prices | -113976.73 | -86.06 | 0.65 | 19.60 | 190.83 |
| Low MP and medium LP | -9921.09 | -39.06 | 21.40 | 41.16 | 173.63 |
| Low MP and high LP | -8513.39 | -40.31 | 38.71 | 62.52 | 787.29 |
| Medium MP and low LP | -1969.68 | -74.81 | -4.63 | 12.34 | 119.20 |
| Medium MP and LP | -5286.69 | -62.31 | 15.39 | 33.04 | 96.95 |
| Medium Mp and high LP | -23942.57 | -60.11 | 32.30 | 54.01 | 94.24 |
| High MP and low LP | -9762.58 | -101.03 | -12.00 | 7.50 | 1267.03 |
| High MP and medium LP | -4617.41 | -58.24 | 8.26 | 27.37 | 689.02 |
| High MP and LP | -4234.33 | -44.08 | 26.61 | 48.78 | 948.13 |
netbenefit_lc_all <- cbind(data.frame(baseline_lc, low_mp, med_mp,high_mp, low_lp, med_lp, high_lp, low_mp_lp,low_mp_med_lp, low_mp_high_lp,
med_mp_low_lp, med_mp_lp, med_mp_high_lp, high_mp_low_lp, high_mp_med_lp, high_mp_lp))
netbenefit_lc_all %>%
pivot_longer(cols = everything(), names_to = "scenarios", values_to = "net_benefit_lc") %>%
group_by(scenarios) %>%
mutate(scenarios = factor(scenarios, levels = c("baseline_lc","low_mp", "med_mp", "high_mp",
"low_lp", "med_lp", "high_lp",
"low_mp_lp", "low_mp_med_lp",
"low_mp_high_lp", "med_mp_low_lp",
"med_mp_lp", "med_mp_high_lp",
"high_mp_low_lp", "high_mp_med_lp",
"high_mp_lp"))) %>%
ggplot(aes(x=scenarios, y= net_benefit_lc))+
stat_summary(geom = 'col', fun = median)+
xlab("Market price scenarios") + ylab("Median net benefit for land consolidation (million KES) ")+
theme_minimal()+ coord_flip()+
# theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
scale_x_discrete(labels =c("low_mp" = "Low MP",
"med_mp" = "Medium MP",
"high_mp" = "High MP",
"low_lp" = "Low LP",
"med_lp" = "Medium LP",
"high_lp" = "High LP",
"low_mp_lp" = "Low MP and LP",
"low_mp_med_lp" = "Low MP and Medium LP",
"low_mp_high_lp" = "Low MP and high LP",
"med_mp_low_lp" = "Medium MP and low LP",
"med_mp_lp" = "Medium MP and LP",
"med_mp_high_lp" = "Medium MP and high LP",
"high_mp_low_lp" = "High MP and low LP",
"high_mp_med_lp" = "High MP and Medium LP",
"high_mp_lp" = "High MP and LP",
"baseline_lc" = "Baseline")) 
Figure 3. The median values of the decision outcome of land consolidation when maize and land prices are varied. These outcomes are based on 1000 model runs of land consolidation model function.
ggsave("figures/fig_3_net_benefit_scenarios.png", width = 9, height = 5) We also calculated the absolute and relative change of the net benefit for land consolidation of different market scenarios.
# mcSimulation_results$y[,"net_benefit_lc_per_ha"] <- mcSimulation_results$y$net_benefit_lc_kes/mcSimulation_results$x$ha_per_hh
# mcSimulation_results$y[,"net_benefit_lc_per_ha_MKes"] <- mcSimulation_results$y$net_benefit_lc_per_ha/1000000
# baseline_lc <- mcSimulation_results$y$net_benefit_lc_per_ha_MKes
netbenefit_lc <- cbind(data.frame(baseline_lc, low_mp, med_mp,high_mp, low_lp, med_lp, high_lp, low_mp_lp,low_mp_med_lp, low_mp_high_lp,
med_mp_low_lp, med_mp_lp, med_mp_high_lp, high_mp_low_lp, high_mp_med_lp, high_mp_lp))
netbenefit_lc[,"low_mp_change"] <- netbenefit_lc$low_mp-netbenefit_lc$baseline_lc
netbenefit_lc[,"med_mp_change"] <- netbenefit_lc$med_mp-netbenefit_lc$baseline_lc
netbenefit_lc[,"high_mp_change"] <-netbenefit_lc$high_mp-netbenefit_lc$baseline_lc
netbenefit_lc[,"low_lp_change"] <-netbenefit_lc$low_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"med_lp_change"]<-netbenefit_lc$med_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"high_lp_change"]<-netbenefit_lc$high_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"low_mp_lp_change"]<-netbenefit_lc$low_mp_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"low_mp_med_lp_change"]<-netbenefit_lc$low_mp_med_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"low_mp_high_lp_change"]<-netbenefit_lc$low_mp_high_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"med_mp_low_lp_change"]<-netbenefit_lc$med_mp_low_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"med_mp_lp_change"]<-netbenefit_lc$med_mp_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"med_mp_high_lp_change"]<-netbenefit_lc$med_mp_high_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"high_mp_low_lp_change"]<-netbenefit_lc$high_mp_low_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"high_mp_med_lp_change"]<-netbenefit_lc$high_mp_med_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,"high_mp_lp_change"]<-netbenefit_lc$high_mp_lp-netbenefit_lc$baseline_lc
netbenefit_lc[,grep("_change$", colnames(netbenefit_lc))] %>%
pivot_longer(cols = everything(), names_to = "market_prices", values_to = "lc_netbenefit_change") %>%
group_by(market_prices) %>%
mutate(market_prices = factor(market_prices, levels = c("low_mp_change", "med_mp_change", "high_mp_change",
"low_lp_change", "med_lp_change", "high_lp_change",
"low_mp_lp_change", "low_mp_med_lp_change",
"low_mp_high_lp_change", "med_mp_low_lp_change",
"med_mp_lp_change", "med_mp_high_lp_change",
"high_mp_low_lp_change", "high_mp_med_lp_change",
"high_mp_lp_change"))) %>%
ggplot(aes(x=market_prices, y = lc_netbenefit_change))+
stat_summary(geom = "col", fun = median) +
xlab("Market price scenarios") + ylab("Median change in net benefit for land consolidation \n (million KES)")+
coord_flip()+ theme_minimal()+
scale_x_discrete(labels =c("low_mp_change" = "Low MP",
"med_mp_change" = "Medium MP",
"high_mp_change" = "High MP",
"low_lp_change" = "Low LP",
"med_lp_change" = "Medium LP",
"high_lp_change" = "High LP",
"low_mp_lp_change" = "Low MP and LP",
"low_mp_med_lp_change" = "Low MP and Medium LP",
"low_mp_high_lp_change" = "Low MP and high LP",
"med_mp_low_lp_change" = "Medium MP and low LP",
"med_mp_lp_change" = "Medium MP and LP",
"med_mp_high_lp_change" = "Medium MP and high LP",
"high_mp_low_lp_change" = "High MP and low LP",
"high_mp_med_lp_change" = "High MP and Medium LP",
"high_mp_lp_change" = "High MP and LP" ))
Figure 4. Change in net benefit for land consolidation given by the difference between outcomes of individual scenarios and baseline. Based on 1000 model runs of the land consolidation model function.
# net_change<-netbenefit_lc[,grep("_change$", colnames(netbenefit_lc))] %>%
# pivot_longer(., cols = everything(), names_to = "market_prices", values_to = "lc_netbenefit_change") %>%
# # mutate(lower_bound = quantile(lc_netbenefit_change, 0.05, na.rm =T),
# # upper_bound = quantile(lc_netbenefit_change, 0.95, na.rm = T)) %>%
# # filter(lc_netbenefit_change>=lower_bound &lc_netbenefit_change<=upper_bound) %>%
# group_by(market_prices) %>%
# summarize(min_quantile = quantile(lc_netbenefit_change, 0.05, na.rm = T),
# median_value = median(lc_netbenefit_change, na.rm = TRUE),
# upp_quantile = quantile(lc_netbenefit_change, 0.95, na.rm = T)) %>%
# pivot_longer(cols=c( "min_quantile", 'median_value',"upp_quantile"),
# names_to = "quantile_levels",values_to = 'change')
#
# ggplot(net_change, aes(x=market_prices, y =change, group = factor(quantile_levels))) +
# geom_point(aes(colour = factor(quantile_levels)))+ geom_path(aes(colour =factor( quantile_levels)))+
# coord_radial(start = 0.25 * pi, end = 1.75 * pi) +
# guides(
# theta = guide_axis_theta(angle = 90),
# r = guide_axis(angle = 0) )+
# scale_x_discrete(labels =c("low_mp_change" = "Low MP",
# "med_mp_change" = "Medium MP",
# "high_mp_change" = "High MP",
# "low_lp_change" = "Low LP",
# "med_lp_change" = "Medium LP",
# "high_lp_change" = "High LP",
# "low_mp_lp_change" = "Low MP and LP",
# "low_mp_med_lp_change" = "Low MP and Medium LP",
# "low_mp_high_lp_change" = "Low MP and high LP",
# "med_mp_low_lp_change" = "Medium MP and low LP",
# "med_mp_lp_change" = "Medium MP and LP",
# "med_mp_high_lp_change" = "Medium MP and high LP",
# "high_mp_low_lp_change" = "High MP and low LP",
# "high_mp_med_lp_change" = "High MP and Medium LP",
# "high_mp_lp_change" = "High MP and LP" )) +
# scale_y_continuous(labels = label_number(suffix = " Million KES")) +
# theme_minimal() +
# labs(x = NULL, y = NULL) +
# scale_colour_manual(name = "",
# labels = c("Median", "5 % quantile","95 % quantile"),
# values = c("deeppink", "limegreen","#1a80bb"))+
# theme(legend.position = c(0.4, 0.97))
ggsave("figures/Fig_4_change_in_decision.png", width = 9, height = 5)netbenefit_lc[,"low_mp_change_relative"] <- (netbenefit_lc$low_mp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"med_mp_change_relative"] <- (netbenefit_lc$med_mp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"high_mp_change_relative"] <-(netbenefit_lc$high_mp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"low_lp_change_relative"] <-(netbenefit_lc$low_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"med_lp_change_relative"]<-(netbenefit_lc$med_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"high_lp_change_relative"]<-(netbenefit_lc$high_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"low_mp_lp_change_relative"]<-(netbenefit_lc$low_mp_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"low_mp_med_lp_change_relative"]<-(netbenefit_lc$low_mp_med_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"low_mp_high_lp_change_relative"]<-(netbenefit_lc$low_mp_high_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"med_mp_low_lp_change_relative"]<-(netbenefit_lc$med_mp_low_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"med_mp_lp_change_relative"]<-(netbenefit_lc$med_mp_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"med_mp_high_lp_change_relative"]<-(netbenefit_lc$med_mp_high_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"high_mp_low_lp_change_relative"]<-(netbenefit_lc$high_mp_low_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"high_mp_med_lp_change_relative"]<-(netbenefit_lc$high_mp_med_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
netbenefit_lc[,"high_mp_lp_change_relative"]<-(netbenefit_lc$high_mp_lp-netbenefit_lc$baseline_lc)/netbenefit_lc$baseline_lc*100
pc_change<-netbenefit_lc[,grep("_change_relative$", colnames(netbenefit_lc))] %>%
pivot_longer(., cols = everything(), names_to = "market_prices", values_to = "lc_netbenefit_pc_change") %>%
group_by(market_prices) %>%
mutate(market_prices = factor(market_prices, levels = c("low_mp_change_relative", "med_mp_change_relative", "high_mp_change_relative",
"low_lp_change_relative", "med_lp_change_relative", "high_lp_change_relative",
"low_mp_lp_change_relative", "low_mp_med_lp_change_relative",
"low_mp_high_lp_change_relative", "med_mp_low_lp_change_relative",
"med_mp_lp_change_relative", "med_mp_high_lp_change_relative",
"high_mp_low_lp_change_relative", "high_mp_med_lp_change_relative",
"high_mp_lp_change_relative")))
ggplot(pc_change, aes(x=market_prices, y = lc_netbenefit_pc_change)) +
stat_summary(geom = "col", fun = median)+ theme_bw()+
scale_x_discrete(labels =c("med_mp_lp_change_relative" = "Medium MP and LP",
"med_mp_low_lp_change_relative" = "Medium MP and low LP",
"med_mp_high_lp_change_relative" = "Medium MP and high LP",
"med_mp_change_relative" = "Medium MP",
"med_lp_change_relative" = "Medium LP",
"low_lp_change_relative" = "Low LP",
"low_mp_med_lp_change_relative" = "Low MP and Medium LP",
"low_mp_lp_change_relative" = "Low MP and LP",
"low_mp_high_lp_change_relative" = "Low MP and high LP",
"low_mp_change_relative" = "Low MP",
"high_mp_med_lp_change_relative" = "High MP and Medium LP",
"high_mp_lp_change_relative" = "High MP and LP",
"high_mp_low_lp_change_relative" = "High MP and low LP",
"high_mp_change_relative" = "High MP",
"high_lp_change_relative" = "High LP"))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
xlab("Market price scenarios") + ylab("Median relative change of net benefit for land consolidation (%)")
Figure S5. Relative change in percent of the net benefit of land consolidation (decision). Based on 1000 moidel runs of the land consolidation model function.
pc_plot <- ggplot(data=pc_change, aes(x=factor(market_prices), y =lc_netbenefit_pc_change))+
geom_violin(draw_quantiles= c(0.5),fill = "turquoise")+
ylab("Relative change (%)") + xlab("Market price scenarios")+
theme_bw() + coord_flip() + ggtitle("Relative change")+
scale_x_discrete(labels =c("low_mp_change_relative" = "Low MP",
"med_mp_change_relative" = "Medium MP",
"high_mp_change_relative" = "High MP",
"low_lp_change_relative" = "Low LP",
"med_lp_change_relative" = "Medium LP",
"high_lp_change_relative" = "High LP",
"low_mp_lp_change_relative" = "Low MP and LP",
"low_mp_med_lp_change_relative" = "Low MP and Medium LP",
"low_mp_high_lp_change_relative" = "Low MP and high LP",
"med_mp_low_lp_change_relative" = "Medium MP and low LP",
"med_mp_lp_change_relative" = "Medium MP and LP",
"med_mp_high_lp_change_relative" = "Medium MP and high LP",
"high_mp_low_lp_change_relative" = "High MP and low LP",
"high_mp_med_lp_change_relative" = "High MP and Medium LP",
"high_mp_lp_change_relative" = "High MP and LP" ))
pc_change2 <- pc_change %>%
group_by(market_prices) %>%
summarize(min_quantile = quantile(lc_netbenefit_pc_change, 0.05, na.rm = T),
median_value = median(lc_netbenefit_pc_change, na.rm = TRUE),
upp_quantile = quantile(lc_netbenefit_pc_change, 0.95, na.rm = T)) %>%
pivot_longer(cols=c( "min_quantile", 'median_value',"upp_quantile"),
names_to = "quantile_levels",values_to = 'pc_change')
ggplot(pc_change2, aes(x=market_prices, y =pc_change, group = factor(quantile_levels))) +
geom_point(aes(colour = factor(quantile_levels)))+ geom_path(aes(colour =factor(quantile_levels)))+
coord_radial(start = 0.25 * pi, end = 1.75 * pi) +
guides(
theta = guide_axis_theta(angle = 90),
r = guide_axis(angle = 0) )+
scale_x_discrete(labels =c("low_mp_change_relative" = "Low MP",
"med_mp_change_relative" = "Medium MP",
"high_mp_change_relative" = "High MP",
"low_lp_change_relative" = "Low LP",
"med_lp_change_relative" = "Medium LP",
"high_lp_change_relative" = "High LP",
"low_mp_lp_change_relative" = "Low MP and LP",
"low_mp_med_lp_change_relative" = "Low MP and Medium LP",
"low_mp_high_lp_change_relative" = "Low MP and high LP",
"med_mp_low_lp_change_relative" = "Medium MP and low LP",
"med_mp_lp_change_relative" = "Medium MP and LP",
"med_mp_high_lp_change_relative" = "Medium MP and high LP",
"high_mp_low_lp_change_relative" = "High MP and low LP",
"high_mp_med_lp_change_relative" = "High MP and Medium LP",
"high_mp_lp_change_relative" = "High MP and LP" )) +
scale_y_continuous(labels = label_number(suffix = "%")) +
theme_minimal() +
labs(x = NULL, y = NULL) +
scale_colour_manual(name = "",
labels = c("Median", "5 % quantile","95 % quantile"),
values = c("deeppink", "limegreen","#1a80bb"))+
theme(legend.position = c(0.4, 0.97))
Figure S5. Relative change in percent of the net benefit of land consolidation (decision). Based on 1000 moidel runs of the land consolidation model function.
8 Unclarified issues
- should the relative change be reported?
References
Luedeling, Eike, Lutz Goehring, Katja Schiffers, Cory Whitney, and
Eduardo Fernandez. 2024. decisionSupport: Quantitative Support of
Decision Making Under Uncertainty. http://www.worldagroforestry.org/.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data
Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.
———. 2023. Tidyverse: Easily Install and Load the Tidyverse. https://tidyverse.tidyverse.org.